Traditionally, two mounds within 100-200m are considered and labelled by archaeologists as a necropolis. Can we determine such clusters computationally? This script shows how to find clusters using mound buffers of 100m and combining this output with Voronoi polygons. Some weakness is in linearly distributed necropoleis, which can escape attention, so clusters need to be composed from bottom up (lowest number of members to the highest number of members)
We have a total of 1090 mounds, extinct and uncertain mound features in Yambol.
# Yambol mounds
features <- readRDS("../output_data/mnds_aggr_27Dec.rds") # dataset is most recent from 27/12/2022
Y_region <- st_read("../data/YamRegion.shp")
## Reading layer `YamRegion' from data source
## `C:\Users\adela\Documents\RStudio\YambolMoundAnalysis\data\YamRegion.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 7 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 431400.8 ymin: 4641684 xmax: 504849.3 ymax: 4727299
## Projected CRS: WGS 84 / UTM zone 35N
head(features)
## Simple feature collection with 6 features and 18 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 441126.7 ymin: 4706447 xmax: 454596 ymax: 4710087
## Projected CRS: WGS 84 / UTM zone 35N
## # A tibble: 6 × 19
## TRAP Source Type LU_Ar…¹ LU_Top Diame…² Heigh…³ Condi…⁴ Princ…⁵
## <dbl> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 9305 <NA> Buri… Scrub Scrub 25 1.8 3 <NA>
## 2 9411 Survey Buri… Scrub Scrub 15 1.5 2 Nodata
## 3 9416 Legacy verification Buri… Annual… Annua… 20 2 3 Agricu…
## 4 9417 Legacy verification Buri… Pasture Pastu… 15 2 4 Looting
## 5 9409 Legacy verification Buri… Annual… Pastu… 30 2 2 Agricu…
## 6 9314 Legacy verification Buri… Pasture Pastu… 30 2 5 Looting
## # … with 10 more variables: geometry <POINT [m]>, distBG [m], distTown [m],
## # elevAster <dbl>, slopeAster <dbl>, aspectAster <dbl>, prom250mbuff <dbl>,
## # TPI <dbl>, TRI <dbl>, rough <dbl>, and abbreviated variable names
## # ¹LU_Around, ²DiameterMax, ³HeightMax, ⁴Condition, ⁵PrincipalSourceOfImpact
features %>% group_by(Type) %>% tally()
## Simple feature collection with 11 features and 2 fields
## Geometry type: GEOMETRY
## Dimension: XY
## Bounding box: xmin: 420515 ymin: 4647347 xmax: 499612.3 ymax: 4722257
## Projected CRS: WGS 84 / UTM zone 35N
## # A tibble: 11 × 3
## Type n geometry
## <chr> <int> <GEOMETRY [m]>
## 1 Burial Mound 1024 MULTIPOINT ((420515 4693580), (420707.1 4690851)…
## 2 Burial Mound? 20 MULTIPOINT ((462092.8 4663311), (465604.9 465350…
## 3 Extinct Burial Mound 188 MULTIPOINT ((420832.5 4689836), (421875.6 469359…
## 4 Extinct Burial Mound? 4 MULTIPOINT ((487113.7 4684048), (488727 4683746)…
## 5 Other 75 MULTIPOINT ((426968.4 4693384), (462171.2 465546…
## 6 Surface Scatter 64 MULTIPOINT ((421045.4 4690265), (424207.4 469496…
## 7 Tell 3 MULTIPOINT ((465676.3 4675875), (467848.2 466240…
## 8 Tell? 1 POINT (483400.7 4697154)
## 9 Uncertain Feature 36 MULTIPOINT ((437475.3 4690926), (445710.4 469931…
## 10 Uncertain Feature? 1 POINT (480327.6 4699085)
## 11 Uncertain Mound 68 MULTIPOINT ((422391.7 4698832), (429156.7 469547…
features <- features %>%
#st_drop_geometry() %>%
dplyr::mutate(Condition = str_extract(Condition, "\\d")) %>%
dplyr::mutate(Condition = case_when(Condition == 0 ~ "NA",
Condition == 6 ~ "5",
Condition != 0 ~ Condition)) %>%
dplyr::mutate(Condition = as.factor(Condition)) %>%
dplyr::mutate(Type= gsub("\\?","",Type)) %>%
dplyr::mutate(HeightMax = as.numeric(HeightMax)) %>%
dplyr::mutate(DiameterMax = as.numeric(DiameterMax))
## Warning in mask$eval_all_mutate(quo): NAs introduced by coercion
mnds <- features %>%
dplyr::filter(Type == "Burial Mound" | Type == "Extinct Burial Mound" | Type == "Uncertain Mound")
mnds
## Simple feature collection with 1304 features and 18 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 420515 ymin: 4647591 xmax: 499612.3 ymax: 4722135
## Projected CRS: WGS 84 / UTM zone 35N
## # A tibble: 1,304 × 19
## TRAP Source Type LU_Ar…¹ LU_Top Diame…² Heigh…³ Condi…⁴ Princ…⁵
## * <dbl> <chr> <chr> <chr> <chr> <dbl> <dbl> <fct> <chr>
## 1 9305 <NA> Buri… Scrub Scrub 25 1.8 3 <NA>
## 2 9411 Survey Buri… Scrub Scrub 15 1.5 2 Nodata
## 3 9416 Legacy verificati… Buri… Annual… Annua… 20 2 3 Agricu…
## 4 9417 Legacy verificati… Buri… Pasture Pastu… 15 2 4 Looting
## 5 9409 Legacy verificati… Buri… Annual… Pastu… 30 2 2 Agricu…
## 6 9314 Legacy verificati… Buri… Pasture Pastu… 30 2 5 Looting
## 7 9315 Legacy verificati… Buri… Pasture Pastu… 20 1 2 Looting
## 8 9415 Legacy verificati… Buri… Annual… Annua… 25 1 5 Agricu…
## 9 9418 Legacy verificati… Buri… Pasture Pastu… 50 6 2 Looting
## 10 9319 Legacy verificati… Buri… Annual… Pastu… 15 1 3 Agricu…
## # … with 1,294 more rows, 10 more variables: geometry <POINT [m]>, distBG [m],
## # distTown [m], elevAster <dbl>, slopeAster <dbl>, aspectAster <dbl>,
## # prom250mbuff <dbl>, TPI <dbl>, TRI <dbl>, rough <dbl>, and abbreviated
## # variable names ¹LU_Around, ²DiameterMax, ³HeightMax, ⁴Condition,
## # ⁵PrincipalSourceOfImpact
mnds_buff200 <- st_buffer(mnds,200)
How many mounds have one or more neighbor(s) within 100 m? Let’s
calculate the number of neighbors via a sdgp matrix with
st_intersects() and then sum() and
filter() help us determine their identity. We can review
their location with mapview()
library(mapview)
mapview(mnds_buff200) + mapview(mnds)
test <- st_intersects(mnds_buff200, mnds_buff200$geometry)
test0 <- lengths(st_intersects(mnds_buff200, mnds_buff200$geometry))>0 # intersects at least one - everything as each intersects itself
test1 <- lengths(st_intersects(mnds_buff200, mnds_buff200$geometry))>1 # 848
test2 <- lengths(st_intersects(mnds_buff200, mnds_buff200$geometry))>3 # 361 mounds in in clusters of 2+ mounds
test5 <- lengths(st_intersects(mnds_buff200, mnds_buff200$geometry))>6 # 147 features 95-97
test6 <- lengths(st_intersects(mnds_buff200, mnds_buff200$geometry))>7 # 119
test7 <- lengths(st_intersects(mnds_buff200, mnds_buff200$geometry))>8 # 87
test10 <- lengths(st_intersects(mnds_buff200, mnds_buff200$geometry))>11 # 30
test11 <- lengths(st_intersects(mnds_buff200, mnds_buff200$geometry))>12 # none
test21 <- lengths(st_intersects(mnds_buff200, mnds_buff200$geometry))>20 # none
# How many clusters of 2 or more mounds are there?
sum(test2)
## [1] 361
# How many clusters of 10?
sum(test7)
## [1] 87
# Which mounds are in clusters of 10?
mnds %>% filter(test10)
## Simple feature collection with 30 features and 18 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 467542.2 ymin: 4683521 xmax: 478932.4 ymax: 4721394
## Projected CRS: WGS 84 / UTM zone 35N
## # A tibble: 30 × 19
## TRAP Source Type LU_Ar…¹ LU_Top Diame…² Heigh…³ Condi…⁴ Princ…⁵
## * <dbl> <chr> <chr> <chr> <chr> <dbl> <dbl> <fct> <chr>
## 1 8022 RS:SITE Buri… Annual… Annua… 20 0.7 5 Agricu…
## 2 8024 RS:SITE Exti… Annual… Annua… 20 0.5 5 Agricu…
## 3 8023 RS:SITE Exti… Annual… Annua… 25 0.1 4 Agricu…
## 4 8025 RS:SITE Exti… Annual… Annua… 20 0.5 5 Agricu…
## 5 8026 RS:FNEG Exti… Annual… Annua… 20 0.3 5 Agricu…
## 6 8027 RS:FNEG Buri… Annual… Annua… 30 1 4 Agricu…
## 7 8028 RS:FNEG Buri… Annual… Annua… 40 1 4 Agricu…
## 8 8021 RS:FNEG Exti… Annual… Annua… 20 0.3 5 Agricu…
## 9 8360 Survey Buri… Annual… Annua… 40 1.4 3 Agricu…
## 10 8365 Legacy verificati… Buri… Annual… Annua… 35 1.1 3 Agricu…
## # … with 20 more rows, 10 more variables: geometry <POINT [m]>, distBG [m],
## # distTown [m], elevAster <dbl>, slopeAster <dbl>, aspectAster <dbl>,
## # prom250mbuff <dbl>, TPI <dbl>, TRI <dbl>, rough <dbl>, and abbreviated
## # variable names ¹LU_Around, ²DiameterMax, ³HeightMax, ⁴Condition,
## # ⁵PrincipalSourceOfImpact
# View the clusters of 5 to 10:
# We add the residual mounds so that we see all the neighbors as not all get tagged!
mnds %>% filter(test5) %>%
mapview() + mapview(mnds, cex = 0.1, zcol = "Type")
At least 10 mounds come out as having 10 neighbors within 100m, however, each of these clusters contains duplicates - mounds which were registered repeatedly under different numbers in successive seasons. Not all of the mounds inside these clusters are marked as members of 10-mound-cluster, as only the central ones have all 10 within their buffer. The mounds at the far edges of the cluster do not get flagged as their buffers do not contain all the cluster fellows.
Furthermore, at least two clusters contain duplicates: - In the cluster W of Karavelovo, the three duplicates are 8022, 8023, 8024, and 8028 which were rudimentarily recorded in 2009 as part of remote sensing ground truthing, and again in 2010 as 9594, 9593, and 9595 and 9592 (from 2010+). IN this case 2009 can be eliminated. This means this cluster has 4 mounds less than counted. - There is another cluster near the quarry west of the Mogila village, with the following five duplicates: 8346, 8347, 8348, 8349, 8350 from 20XX to 9222, 9223, 9224, 9225, 9226. These need to be reviewed.
# Sanity check on the duplicates in the other clusters
mnds[95,]
## Simple feature collection with 1 feature and 18 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 464359.3 ymin: 4706914 xmax: 464359.3 ymax: 4706914
## Projected CRS: WGS 84 / UTM zone 35N
## # A tibble: 1 × 19
## TRAP Source Type LU_Ar…¹ LU_Top Diame…² Heigh…³ Condi…⁴ Princ…⁵
## <dbl> <chr> <chr> <chr> <chr> <dbl> <dbl> <fct> <chr>
## 1 9222 Legacy verification Buri… Pasture Pastu… 20 2 3 Looting
## # … with 10 more variables: geometry <POINT [m]>, distBG [m], distTown [m],
## # elevAster <dbl>, slopeAster <dbl>, aspectAster <dbl>, prom250mbuff <dbl>,
## # TPI <dbl>, TRI <dbl>, rough <dbl>, and abbreviated variable names
## # ¹LU_Around, ²DiameterMax, ³HeightMax, ⁴Condition, ⁵PrincipalSourceOfImpact
test[[95]] # buffer around point 95 has allegedly over 6 neighbours
## [1] 95 96 97 347 348 609 610 613 614 616 617
#[1]95 96 97 347 348 610 613 614 616 617
mapview(mnds[c(95:97,347:348,613:614,616:617),])
# Duplicates - Mogila cluster!
dplyr::distinct(mnds)
## Simple feature collection with 1304 features and 18 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 420515 ymin: 4647591 xmax: 499612.3 ymax: 4722135
## Projected CRS: WGS 84 / UTM zone 35N
## # A tibble: 1,304 × 19
## TRAP Source Type LU_Ar…¹ LU_Top Diame…² Heigh…³ Condi…⁴ Princ…⁵
## <dbl> <chr> <chr> <chr> <chr> <dbl> <dbl> <fct> <chr>
## 1 9305 <NA> Buri… Scrub Scrub 25 1.8 3 <NA>
## 2 9411 Survey Buri… Scrub Scrub 15 1.5 2 Nodata
## 3 9416 Legacy verificati… Buri… Annual… Annua… 20 2 3 Agricu…
## 4 9417 Legacy verificati… Buri… Pasture Pastu… 15 2 4 Looting
## 5 9409 Legacy verificati… Buri… Annual… Pastu… 30 2 2 Agricu…
## 6 9314 Legacy verificati… Buri… Pasture Pastu… 30 2 5 Looting
## 7 9315 Legacy verificati… Buri… Pasture Pastu… 20 1 2 Looting
## 8 9415 Legacy verificati… Buri… Annual… Annua… 25 1 5 Agricu…
## 9 9418 Legacy verificati… Buri… Pasture Pastu… 50 6 2 Looting
## 10 9319 Legacy verificati… Buri… Annual… Pastu… 15 1 3 Agricu…
## # … with 1,294 more rows, 10 more variables: geometry <POINT [m]>, distBG [m],
## # distTown [m], elevAster <dbl>, slopeAster <dbl>, aspectAster <dbl>,
## # prom250mbuff <dbl>, TPI <dbl>, TRI <dbl>, rough <dbl>, and abbreviated
## # variable names ¹LU_Around, ²DiameterMax, ³HeightMax, ⁴Condition,
## # ⁵PrincipalSourceOfImpact
duplicates <- lengths(st_equals(mnds))>1 # produces zero duplicates, as they do not have identical locations, but are withing 10-20m of one another. A spatial filter on close proximity is needed
# Close mounds
mnds[c(95:97,347:348,613:614,616:617),]
## Simple feature collection with 9 features and 18 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 464228.7 ymin: 4706914 xmax: 464404.7 ymax: 4706973
## Projected CRS: WGS 84 / UTM zone 35N
## # A tibble: 9 × 19
## TRAP Source Type LU_Ar…¹ LU_Top Diame…² Heigh…³ Condi…⁴ Princ…⁵
## <dbl> <chr> <chr> <chr> <chr> <dbl> <dbl> <fct> <chr>
## 1 9222 Legacy verification Buri… Pasture Pastu… 20 2 3 Looting
## 2 9225 Legacy verification Buri… Pasture Pastu… 15 1 2 None
## 3 9226 Legacy verification Buri… Pasture Pastu… 20 1.5 3 Looting
## 4 9223 Survey Buri… Pasture Pastu… 15 1.5 4 Looting
## 5 9224 Survey Buri… Pasture Pastu… 20 1.5 3 None
## 6 8347 Legacy verification Buri… Pasture Pastu… 17 1.2 1 Post-d…
## 7 8348 Legacy verification Buri… Pasture Pastu… 15 1.5 2 Looting
## 8 8349 Legacy verification Buri… Pasture Pastu… 17 1.75 2 Looting
## 9 8350 Legacy verification Buri… Pasture Pastu… 18 1.2 2 Post-d…
## # … with 10 more variables: geometry <POINT [m]>, distBG [m], distTown [m],
## # elevAster <dbl>, slopeAster <dbl>, aspectAster <dbl>, prom250mbuff <dbl>,
## # TPI <dbl>, TRI <dbl>, rough <dbl>, and abbreviated variable names
## # ¹LU_Around, ²DiameterMax, ³HeightMax, ⁴Condition, ⁵PrincipalSourceOfImpact
There is ast_voronoi() function but when using it on
points, one must st_combine() the points first (as with
convex hull) so as to create a mesh covering all the points.
# Create Voronyi geometries
vor <- st_voronoi(st_combine(mnds))
# st_voronoi returns a GEOMETRYCOLLECTION,
# some plotting methods can't use a GEOMETRYCOLLECTION.
# this returns polygons instead
mnd_vor_poly <- st_collection_extract(vor)
# Crop the polygons to Yambol region boundary
vor <- st_intersection(Y_region, mnd_vor_poly)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
# Visualize
mapview(vor)+ mapview(mnds %>% filter(test5)) + mapview(mnds, cex = 0.1)
Nice overview of where some of the medium sized clusters are. Next step: mark the clusters by 1-10 depending on how many neighbors each mound has and produce a colour-coded version so one can see the clusters better.